Quantum regression theorem for non-Markovian Lindblad equations 
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We find the conditions under which a quantum regression theorem can be assumed valid for 
non-Markovian master equations consisting in Lindblad superoperators with memory kernels. Our 
considerations are based on a generalized Born-Markov approximation, which allows us to obtain 
our results from an underlying Hamiltonian description. We demonstrate that a non-Markovian 
quantum regression theorem can only be granted in a stationary regime if the dynamics satisfies a 
quantum detailed balance condition. As an example we study the correlations of a two level system 
embedded in a complex structured reservoir and driven by an external coherent field. 
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I. INTRODUCTION 

In many areas of physics one is confronted with the de- 
scription of small quantum systems interacting with an 
uncontrollable environment. This situation is well under- 
stood when the reduced system dynamics follows a (com- 
pletely positive) Markovian evolution [E II II S Li l| • 

One of the cornerstone of the theory of Markovian 
open quantum systems is the quantum regression theo- 
rem (QRT). This theorem, originally proposed by Lax 0, 
allows to calculate multiple-time operators correlation 
functions from the knowledge of single-time expectation 
values, which in turn implies the knowledge of the density 
matrix evolution |J, li| ■ The importance of this the- 
orem comes from the physical information contained in 
the operator correlations. In fact, in a stationary regime, 
it is possible to relate the Fourier transform of these ob- 
jects with the spectrum of the decay process Q ■ Further- 
more, in radiant systems, the statistic of the scattered 
field can be described through system operator correla- 
tions HSU. 

Another central cornerstone of non-equilibrium quan- 
tum Markovian dynamics is the quantum detailed bal- 
ance condition, which imposes severe symmetry prop- 
erties on the operator correlations structure. While in 
classical stochastic processes this condition has a clear 
meaning in terms of transitions between the available 
states of the system gj, |9|| , in quantum dissipative sys- 
tems this condition relies in the time reversal property of 
the underl ying stationary microscopic Hamiltonian evo- 
lution [13, HI El El El El- The breakdown of this 
condition has direct experimental implications |l6j| . 

Although the applicability of the Markovian approxi- 
mation range over many physical situations 0, H H 0, 
H 0, there exist several real systems whose dynamics 
present strong departures from it. Remarkable examples 
are anomalous intermittent fluorescence in quantum dots 



[171 Il8l Il9t 1 2 ("| , the presence of 1/ f noise in phase and 
charge su p er conducting qubits [2l[, [23, and band gap 
materials 1H 

Consistently with the existence of experimental situa- 
tions that can not be described by a Markovian evolu- 
tion, in the context of different approaches recent effort 
was dedicated to characterize non-Markovian operator 
correlation dynamics [25L l2rj |27| . 

While the description of non-Markovian processes may 
depends on each specific situation, there exists an in- 
creasing interest in describing these kind of processes by 
introducin g m emory contributions in standard Lindblad 
evolutions [28|, |2JJ, yCJ, IU |H IH yj, |35|, |3(| . This proce- 
dure provide easy manageable equations. Nevertheless, 
this technique does not have associated a rule for calcu- 
lating operator correlations. 

In this paper we explore the possibility of establishing 
a QRT for non-Markovian master equations that can be 
cast in the form of Lindblad equations with memory con- 
tributions m, m m m m m m m . we wm base 

our considerations in a generalized Born-Markov approx- 
imation (GBMA) j3||, which allows us to develop our 
results from an underlying microscopic Hamiltonian de- 
scription. 

The paper is outlined as follows. In Section II we 
review the derivation of the GBMA from a full micro- 
scopic description. Based on this approach, in Section III 
we search the conditions under which a non-Markovian 
QRT can be established. In Section IV we relate the 
non-Markovian QRT with a detailed balance condition. 
In Section V we exemplify our theoretical results by an- 
alyzing the correlation dynamics of a two level system 
embedded in a complex thermal environment described 
in a GBMA. In section VI we give the conclusions. 



II. 



GENERALIZED BORN-MARKOV 
APPROXIMATION 



* present address 



The GBMA applies for complex structured environ- 
ments whose action over the system can be well approx- 
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imated by a direct sum of sub-reservoirs, each one being 
able to induce by itself a Markovian system dynamics. 
Under this condition, the system evolution can be writ- 
ten as a Lindblad equation characterized by a random 
dissipative rate [3t| . 

Here, we review the microscopic derivation of this ap- 
proximation by using a well known projector operator 
technique |37l l38| . This equivalent derivation is useful for 
clarifying that the GBMA is not restricted to a second 
order approximation. In fact, the projector technique 
provides a rigorous procedure that allows to obtain the 
system dynamics up to any desired order in the interac- 
tion Hamiltonian. 

We assume a full microscopic Hamiltonian description 
of the interaction of a system S with its environment B 



Ht — Hs + Hb + Hj. 



(1) 



Here, H$ and Hb correspond to the system and bath 
Hamiltonians respectively. The term Hj = q$ <8> Qb de- 
scribes their mutual interaction, with the operators qs 
and Qb acting on the system and bath Hilbert spaces 
respectively. 

The system density matrix follows after eliminating 
the environment degrees of freedom, ps(t) = Tr B{pr(t)}, 
where the total density matrix prit) evolves as 

d f^L = ^[H T ,pr(t)] = C T [ PT {t)]. (2) 

The GBMA j3|| can be derived by introducing the pro- 
jector V defined by 



Tp T (t) = J2pR(t)®Z R , 



where ^r is given by 



^R 



HrPbHr, 



(3) 



(4) 



with pb being the stationary state of the bath, while the 
system states pn{t) are defined by 



PR® 



TTB{Tl R p T (t)n R } 

Tr B {n_ R/ 9 B n_ R } 



(5) 



Here, we have introduced a set of projectors = 
S{ £R } \ e R) ( e r \ i which provides an orthogonal decompo- 
sition of the unit operator [Ig] in the Hilbert space of 
the bath, J^r^r = I B , with UrUr> = Ur5r, w . The 
full set of states \€r) corresponds to the base where ps 
is diagonal, which implies J^r^R — Pb- 

It is easy to realize that V 2 = P. In physical terms, 
this projector takes in account that each bath-subspace 
associated to the projectors Hr induces a different sys- 
tem dynamics, each one represented by the states PR(t). 
On the other hand, notice that the standard projector 
Vpr(t) = Tr B {pT(t)} ® pb = ps{t) ® Pb WK E3, is 
recuperated when all the states PR(t) have the same dy- 
namics. 



From Eq. ©, the system density matrix follows as 
ps(t) = Tr B {Vp T (t)} = Y, P RPR(t) = (M*)>- ( 6 ) 



This equation defines the system state as an average over 
the density matrixes PR(t), each one participating with 
weight Pr. These parameters arc defined by the weight 
of each subspace in the full stationary bath state 

Pr = Tr B {S R } = Tt b {^-rPb} = ^ ( 6R \ <° s l £fl ) ' ( 7 ) 

which in consequence satisfy ^2rPr. = 1. 

By writing the evolution Eq. J2J) in an interaction rep- 
resentation, and splitting the full dynamics in contribu- 
tions Ppx(t) and Qpr(t), where Q = 1 — V, up to second 
order in the interaction Hamiltonian it follows 37, 38] 



dP PT {t) 
dt 



dt'PCT(t)CT(t')Pp T (t') 



(8) 



where Cx(t) is the total Liouville operator in a interac- 
tion representation. Here, we have assumed an uncorre- 
cted initial state, pr(0) = Ps(0) (g> pb- 

By assuming an interaction Hamiltonian with a direct 
sum structure 



Hr = H r 



Ht 



(9) 



where each term satisfies Hj R = UrHiIIr, from Eq. © 
it follows that each state PR(t), in a Schrodinger repre- 
sentation, evolves as 



dpn{t) 



■{H Sl p R {t)} 



dt' 



1 

dt h L ° ' rr<A " \h; ./ii 

Tr BR {[H lR , [H lR (-t r ),p R (t) ® pb r }}. 



(10) 



with p Blt = Er/Pr, and where Tt Br {»} = Tr B {IL R • 
Hr}. The corresponding initial condition reads pr{0) = 
ps(0), which follows from Eq. (JSJ. Furthermore, in this 
evolution we have introduced a Markov approximation, 
which applies when each bath subspace corresponding to 
the projectors Hr defines a Markovian sub-environment. 

The evolution Eq. I|1(J|I , disregarding transients of the 
order of the sub-bath correlation time, can be always well 
approximated by a Lindblad equation 



dpR.{t) 
dt 



C H { P R{t)]+ 1 RC[pR{t% (11) 



where £#[•] = —(i/H)[Hs, •], and the dissipative contri- 
bution is defined by a Lindblad superoperator Q] 

£[•] = ~ E a <*p{[V a , + [V a ; Vj]). (12) 

a/3 

Here, the set of system operators {V a } and the dimen- 
sionless Hermitian matrix a a p depend on the underly- 
ing microscopic interaction. The rates jr follow from a 
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Fermi golden rule when applied to the manifold of states 
{|efl}} that define each Markovian sub-reservoir. 

While the density matrixes pn{t) follow a Markovian 
evolution, the system state ps(i) evolves with a com- 
pletely positive |l|,[2j non-Markovian evolution, property 
inherited from the random Lindblad structure Eq. Q11JI. 
The average of this equation over the set {7^, Pr} can 
be performed in a Laplace domain, from where it follows 

= C H [ Ps {t)] + J* drh(t r)[p s (r)}, (13) 

where the superoperator L(t) is defined by the relation 



expectation values 



(G R (u) lR £)[.} = (G R (u))L(u)[* 



(14) 



Here, u is the Laplace variable and Gr(u) is the Marko- 
vian propagator of each state pn{t), i.e., Gr(u) = [u — 
(Ch + 7r£)] _1 . Depending on the set {7^, Pr}, which 
specify the complex environment, Eq. (|13|) may lead to 
a reach variety of system decay behaviors as well as to 
many different structures of non-local Lindblad equations. 
In fact, in general L(u) consists in a sum of Lindblad 
terms, each one characterized a different memory kernel. 

The structure of the superoperator L(u ) can be sim- 
plified in an effective approximation [36J, which con- 
sists of discarding the dependence introduced by the 
Lindblad superoperator C in the propagator Gr(u), i.e., 
Cr — > —I. From Eq. I|14|) it follows the approximated so- 
lution L(u) ~ K(u — £h)£, which implies the evolution 

~ C H [p s (t)} + J drK(t - r)e^- T)c « C[ Ps {t)]. 

(15) 

with K{u) = (jr{u + Jr)- 1 )((u + T^rV 1 - H the 
time scale of the unitary dynamics is larger than the 
time scale of the memory kernel, the unitary contribu- 
tion can be discarded leading to a single memory Lind- 
blad equation. We remark that structures similar to 
Eq. (I15|l were obtained in the context of other approaches 
IH El |H |H HU. The GBMA, here defined 
through the projector Eq. ©, allows us to associate an 
underlying well defined microscopic description to these 
kind of equations. 



III. QUANTUM REGRESSION THEOREM 

For Markovian master equations the QRT 0,|5, 6] pro- 
vides a direct relation between the evolution of the ex- 
pectation values of system observables and their corre- 
sponding correlation functions. Here we will explore the 
possibility of formulating an equivalent relation when the 
system dynamics can be described through the GBMA. 



Random rate formulation for operators 
correlations 



Let us introduce a complete set of operators {^4^} of 
the system, collected into a vector A, and consider the 



A(t)=Tr SB [A(t)p T (0)], 
as well as the correlation functions 



(16) 



0(t)A(t + t) = Ti SB [0{t)A(t + T)p T (0)], (17) 

where 0(t) is an arbitrary system operator. The time 
dependence of the operators refers to a Heisenberg rep- 
resentation with respect to the total Hamiltonian Eq. JJ), 
i.e., 0(t) = exp[(i/h)tH T ]0(Q)exp[-(i/H)tH T ]. 

From Eq. JHJ, we can write the expectation values as 
an average over the solutions corresponding to each rate 



A(t) = <Tr s [A(0)pfl(*)]> ee (A(t) R ). 



(18) 



In order to work out the operator correlations, 
we first express the total initial density matrix as 
p T (0) = exp[(i/h)tH T ]p T (t)exp[-(i/h)tH T ]. Then, by 
using the cyclic property of the trace, from Eq. I|17fl we 
obtain 



Q{t)A(t + t)= Trs{A(0)Tr s [OsB(r)]}, (19) 



where the operator Osb(t) satisfies 

^O sb (t) = ~[H T) Osb{t)], 



(20) 



with Osb( t )\ t= q — PT(t)O(0). This system-bath opera- 
tor evolves as the total density matrix, Eq. (J2J). On the 
other hand, Eq. (J3J) allows us to write the initial condition 
as O sb {t)\ t=q » J2 R {pR(t)O{0)]®ZR [U. Therefore, 
the reduced dynamics of Osb(t) can also be described 
in a GBMA, which deliver 

Tt b [O sb (t)} = (cxp[(£ H + C R )T}pR(t))O(0), (21) 

where, for shortening the notation we defined Cr = ^rC. 
From Eq. ljl§|) , it follows 



0(t)A(t + r) = (TT S {A(0)e^ +c ^ T [ P R(t)O(0)}}) 

= (0(t)A(t + r) R ). (22) 

This expression is an average over the random set 
{lR, Pr} °f the corresponding Markovian correlation ex- 
pressions 0. This characteristic provides us a central 
result, which allows us to extend the averaging proce- 
dure [Eq. ©] corresponding to the GBMA for operator 
correlations as well. In fact, higher correlations operators 
can also be obtained as an average, over the random rate 
set, of the Markovian expressions corresponding to each 
state PR(t). For example, using the same steps as before, 
for arbitrary system operators 0\ and O2, it is possible 
to obtain 

Oi(i)A(i + r)0 2 (*) - {TY S {Ae^" +c ^ r [0 2 pR(t)0 1 }}} 
= (0^t)A(t~T7)dW) R ), (23) 

which also correspond to an average over the correspond- 
ing Markovian dynamics Q . 
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B. Expectation and correlation evolution 

From the previous result, we can write the evolution of 
both, expectation values and correlations, as an average 
over the random rate set 



-A(t) = (M R A(t) R ), 



(24a) 



— 0(t)A(t + r) = (M R 0(t)A(t + T) R ). (24b) 

Here, the matrix M# acts on the indices of A and is 
defined by the condition 

Tr s {A(£ H + C R )[0}} = M R Tr s {AO}. (25) 

When 7^ is fixed, the evolution equations l|24a|) for ex- 
pectation values and (|24b(l for correlation functions are 
identical, which recovers the QRT for Markovian dynam- 
ics. In the non-Markovian case, however, both equations 
still involve the average over the dissipation rate. 

As for the density matrix (3(| , the averaged evolutions 
can be worked out in the Laplace domain. The expec- 
tation value can be expressed as A(u) — {G R (u)A(0)), 
with the matrix propagator G R (u) = (u + M^) _1 . Af- 
ter introducing the identity operator in the form A(u) = 
(G R (u)(u + Mi?)) -1 (Gr(u)A(O)) , we arrive to the de- 
terministic closed evolution 



-Mt) 



dt j Q 



dt'M(t - t')A(t'). (26a) 



Using a similar procedure, for the correlation we get 
d 



. 0{t)A{t + r) = - / dt 1 M(r - t')0(t)A{t + t') 
dT Jo 



-I(t,r) 



The deterministic kernel matrix 



(26b) 

fulfills the equation 



M(u) =(G R ( U ))- 1 (G R ( U )M R ) ) 
while the inhomogeneous term I(t, t) is defined by 



(27) 



ItMlHGflWl-^GflHOWAff),) - (0(t)A(t) R ). 

(28) 

Besides that Eq. (|24b(l has the same structure as 
Eq. (|24a|l . the inhomogeneous term is only present in 
the correlation evolution, Eq. I)26b|) . I(t, t) arise because 
the initial condition of each contribution in Eq. J24bjl is 
correlated with respect to its propagator. In fact, notice 
that both G R (u) and 0(t)A(t) R depend on j R , which 
implies that theses objects are correlated with respect t o 
the random rate statistics. The dependence of 0{t)A(t) R 
on lR follows from 0{t)A{t) R = Tt s {O(0)A(0)p R (t)}. 
On the other hand, as Eq. (|24a|l is defined with initial 



not depends on j R , which in turn implies that the inho- 
mogeneous term is not present in the averaged evolution 
Eq. (Hi}) . 

Due to the inhomogeneous term I(t, r), the QRT is not 
fulfilled in general. A non-Markovian QRT is only valid 
when this term vanish, which leads to the condition 



(G R (u)0(t)A(t) R ) Q ^ T (G R (u))(0(t)A(t) R ) (29) 

This equality is always satisfied for Markovian dynam- 
ics because the average over the dissipation rate is ab- 
sent. We also note that a non-Markovian QRT can be 
asymptotically valid if the stationary state = p R (oo) 
does not depend on ^y R 40]. In fact, in this situation 
lim^oo 0(t)A(t) R = Tr s {O(0)A(0)p^] is independent 
of 7i?, and then the condition Eq. (|29J) is automatically 
satisfied. However, if the asymptotic state p R depends 
on 7# the inhomogeneous term will contribute at all 
times, even in the asymptotic regime, and the QRT is 
invalidated. The same condition is valid for higher oper- 
ators correlations. 



C. Non-Markovian dynamics 



The evolution Eq. (|26a|l and (|26 b|) can be formally inte- 
grated in the Laplace domain. For the expectation values 
we get 



A(t) =G(i)A(0), 
while for the correlations it follows 



(30a) 



0(t)A(t +t) = G(t)0(*) A(t) + F(t, r). (30b) 
The non-Markovian propagator is defined by 

1 



G(«) 



u + M(u) 

and the extra inhomogeneous term is 



(31) 



F(t,r) = (G R (r)0(t)A(t) R ) - (G R (r))(0(t)A(t) R ). 

(32) 

These expressions explicitly show that the departure 
from condition Eq. I|29|l measures the size of the dynam- 
ical effects which can not be captured by assuming valid 
the QRT. In fact, the QRT is fulfilled only when F(t,r) 
vanishes. 

Eq. I|30al) and l|30b|l are consistent with the averaging 
procedure over Markovian solutions. In fact, they can be 
expressed as A(t) = (G_r(t))A(0), and for the correla- 
tions as 



conditions fixed at t = 0, its initial condition A(0) R does 



0(t)A(t + t) = (G R (r)0(t)A(t) R ), (33) 
which in fact are an average over Markovian solutions. 
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D. Fluctuation operators 

Of special interest is to study the correlation dynamics 
of fluctuation operators, which are defined as the depar- 
ture from expectation values 



SA(t) = A(i)-A(i), 
SO(t) = 0(t)-0(Tj. 



(34a) 
(34b) 



For Markovian dynamics the correlation of these opera- 
tors also satisfies a QRT. These objects are relevant to 
split the spectrum, defined as the Fourier transform of 
the stationary correlations, in a coherent and incoherent 
components [j]. 

From Eq. (|33|) we can write 



0(t)A(t 



(G R (T)dO(t)SA(t + T ) B ) 

+ {W) R Mt + r) R ). (35) 



For Markovian evolutions, in the asymptotic time regime 
(t — > oo), the first contribution can be associated with the 
incoherent spectrum component while the second one, 
after taking the extra limit r — > oo, with the coherent 
spectrum part. After averaging over the random rate, 
these associations remains valid for the non-Markovian 
case. In particular, we note that the coherent component 



lim 0(t)A(t + t) = (0(oo) R A(oo) R ) (36) 



is an average of the corresponding Markovian contribu- 
tions. 

The correlation Eq. (|33(1 can also be written as 



0(t)A(t + r) = SO(t)SA(t + t) + 0(t) A(t + r). (37) 

As for the Markovian case, this expression follows imme- 
diately from the microscopic definition Eq. 117fl . From 
this relation and Eq. Ij35|l . we get 



SO(t)5A(t + t)= G{r)SO(t)SA(t) + 5F(t, r). (38) 

Here, the first contribution follows from the QRT when 
assumed valid for fluctuations operators, and the second 
one measures the departure from it, being defined by 



SF(t,r) = (G R (r)SO(t)dA(t) R )-(G R (T))(SO(t)SA(t) R ) 
+ (0(t) R A(t + T) R ) - (0{t) R )(Mt + r) R ). 

As for operators, in the asymptotic regime the QRT is 
also valid for fluctuation operators if the stationary sate 
p R does not depend on the random rate, which is fact 
implies <5F(oo,t) = 0. 

From Eq. I|37|l and Eq. I]38p. we notice that by assum- 
ing valid the QRT, the coherent spectrum component 
reads 



■ QRT 



lim 0{t)A(t + t) = (0(oo) R )(A(oo)„). (39) 



This expression and Eq. I|36|l indicate that the predictions 
of the QRT will differ from the exact dynamics not only 
in the transient dynamical behaviors but in general also 
in the asymptotic correlation values. In the next sections, 
we will use the difference between these two expressions 
as a measure of the deviation from the validity of the 
QRT in the stationary regime. 



IV. DETAILED BALANCE CONDITION 

In the context of the GBMA, in the previous section 
we have demonstrated that the QRT can be assumed 
valid in an asymptotic regime if the stationary state 
corresponding to each Markovian contribution does not 
depends on the random rate. Here, we will find an equiv- 
alent condition which does not depends on the approx- 
imations used to arrive to the non-Markovian Lindblad 
equation. We will demonstrate that the previous result 
can be associated with a quantum detailed balance condi- 
tion [13, E3, E3, EE E3, US , which in turn is related with 
the microrreversibility of the underlying microscopic dy- 
namics llfl. 



A. Classical conditions 

The concept of detailed balance is well established for 
classical population master equations 



dp n (t) 
dt 



Id 



| ^ gnmPrn (*) - ^ 9mnPn (*) } , (40) 



where "/ c i9nrn define the hopping rates. The classical de- 
tailed balance condition reads 



(oo) (oo), 



(41) 



t — >oo 
r — >oo 



which has an immediate interpretation in terms of the 
available stationary transitions. We note that these rela- 
tions does not depend on the global rate 7 C ;. Thus, they 
impose strong relations between the dimensionless hop- 
ping coefficients {g n m} and the stationary populations 
{p n (oo)}. In particular, it is possible to prove that when 
the stationary state depends on an arbitrary continuous 
parameter e, {p m (oo,e)}, the hopping coefficients must 
also to depend on that parameter, {g mn (£)}- If this is not 
the case, the detailed balance condition is violated [4l| . 
This result can be extended to quantum master equa- 
tions, e being the random rate, establishing a strong re- 
lation between the validity of the QRT for non-Markovian 
dynamics and the detailed balance condition. 



B. Quantum Markovian conditions 

The detailed balance condition can be generalized for 
quantum dynamics from the time reversal property [RH 
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of the stationary system-bath dynamics. For an open 
Markovian system, it can be written as an statement of 
time sym metry for stationary two-time operator correla- 
tions [SIHIlllIl 

lim 0(t + T)A(t) n = lim A(t + r)0(t) H , (42) 

where 0(t) and A(t) represent time-reversed operators 
f4^| . From this equation [l^, it is possible to write an 
equivalent formulation in the Laplace domain as |4flj 

Pr (f l r# A *]= A pr'}- (43) 

We have introduced dual and time reversed superoper- 
ators 0, 0, 0, ^| , which respectively are defined by 

Tr s {OC[p}} = Tr s {pC*[0}}, and by C\0] = C{0] 0. 
Equation 143(1 is equivalent to the conditions 

Pr = Pr, (44a) 
H SP % = p^H Sl (44b) 

7«p££ # M = 7ii«.]. (44c) 

From the second equation and the stationary condition, 
{C H + 1r£}[Pr] = 0, it follows C[p%] = 0. This con- 
dition cannot be satisfied consistently if the stationary 
state p'r depends on jr. In fact, the superoperator C 
does not has a continuous parametrized degenerate null 
eigen-operator. An equivalent conclusion can be obtained 
from the third relation. Then, we deduce that whenever 
Pft depends on the random rate jr the detailed bal- 
ance condition is violated. Therefore, we can affirm that 
if the underlying Markovian evolution of Pij(t) satisfies 
the quantum detailed balance conditions, Eq. \44\I j the 
non-Markovian QRT is valid in the asymptotic regime. 
Equivalently, this statement indicates that when the non- 
Markovian QRT is not fulfilled, the detailed balance con- 
ditions Eq. (|44|l are also not satisfied. 



which must to be valid for any value of the Laplace vari- 
able u m| . We notice that a similar structure also arises 
when formulating the detailed balance condition for non- 
Markovian classical Fokker-Planck equations . 

In contrast to the previous conditions [Eq. (03J], 
Eq. H47JI do not depends on the approximations or for- 
malism used to derive the non-Markovian system dynam- 
ics. In fact, it only depends on the superoperator L(m) 
that defines the density matrix evolution, Eq. (13} • In 
this way we establish a general relation between the non- 
Markovian QRT and the non-Markovian quantum de- 
tailed balance condition. We can affirm that, whenever 
the non-Markovian quantum detailed balance conditions 
Eq. \4 7| ) are satisfied, the non-Markovian QRT is fulfilled 
in the asymptotic regime. The superoperators L# (u) and 
L(it) follow from L(w) after replacing all involved super- 
operators by their dual and time reversed expressions, 
respectively. In the context of the GBMA, they satisfy 
Eq. (|14|l after replacing Cr and C by their dual and time 
reversed expressions. 

As we will exemplify in the next section, a typical sit- 
uation where the non-Markovian QRT is broken, even in 
the stationary regime, is in systems at thermal equilib- 
rium subject to an external perturbation |9j]. In fact, it is 
possible to prove that conditions Eqs. I|47l) are not satis- 
fied when a dissipative dynamics that by itself fulfill the 
detailed balance condition is subject to the action of an 
external Hamiltonian field that does not commutate with 
the system Hamiltonian. Equivalently, in the context of 
the GBMA, the presence of the external perturbation im- 
plies that p^ depends on the random rate, which broke 
the fulfillment of the Markovian conditions Eq. (|44|) . 



V. DECAY IN A STRUCTURED THERMAL 
RESERVOIR 



C. Quantum non-Markovian conditions 

The microrreversibility condition Eq. I|42l) can be triv- 
ially extended to the non-Markovian dynamics as 



lim 0(t + r)A(t) = lim A(t + r)O(i). (45) 

t — >oo t — >oo 

After applying the averaging procedure, from Eqs. 142(1 
and (33Ji wc get the equivalent condition 



Pr 



u -(£*+£#) 




1 



u - (C H + Cr) 



[Pr 



( 46 ) 

When the stationary state does not depends on the ran- 
dom rate, pj? = P's ^1- Q46fl leads to the conditions 



~oo _ 00 

Ps P s ' 



(47a) 

pf{C*+L*(u)}[.] = {£ B + Uu)}\pf*],(m) 



Here we will exemplify our theoretical results by study- 
ing a two level system embedded in a complex structured 
thermal reservoir whose action can be described through 
the GBMA. The system Hamiltonian is 



Ha 



hjjA 



o z + H ext (t), 



(48) 



where Hlua is the difference of energy between the two lev- 
els, denoted by |±), and o~ z is the z-Pauli matrix. H ext (t) 
represent an external time dependent field. 

The dissipative system dynamics can be defined 
through the evolution of the states PR(t), which reads 



= C H [pR{t)] +1 ' R C th [pR(t)] + 2^-u[p R (t)}}. 

(49) 

with £#[•] = —(i/H)[Hs,*\. The influence of the struc- 
tured thermal reservoir is introduced by the Lindblad 
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superoperator 

+^([^ t ,H + k t -^]), (50) 

and an arbitrary set {7 ' r, Pr} of random rates and 
weights, cr^ and a are the raising and lowering opera- 
tors acting on the states |±). The dimensionless con- 
stant n t h defines the temperature T of the environment 
as exp[— huiA/kT] — nth/inth + 1): where k is the Boltz- 
mann constant. We have also considered the action of an 
extra dispersive environment which is introduced by the 
Lindblad superoperator 

£»[•] = + [a z ;a z })/2, (51) 

and the single non-random rate 7$. 

A. Free decay dynamics 

First we analyze the case without the external excita- 
tion, i.e., H ext (t) = 0. 

Density matrix evolution: The evolution of the sys- 
tem density matrix follows from Eq. I|13f) and (|14fl . By 
denoting the matrix elements as 



m _ ( n+(i) *+(<) \ 
w(() -(*.(t) n_(t) J' 



(52) 



in an interaction representation with respect to Hloao~ z /2, 
for the populations we get the evolution 

jU ± (t) = J drJT(t - T){ T n e m + (r) ± n e /n_(r)}, 

(53) 



while for the coherences we obtain 

ft 



d r 

—q> ± (t) = -j drK^{t - r)$ z 



(r). 



The memory kernel functions are defined by 



K{u) 



IB 



1 



" + 7-R/ \U + lR 
1 





7* 



u + 7 fi 



(54) 

(55) 
(56) 



For shortening the notation, we introduced the rates 
7A = 7^(1 + 2n t h) and 7^ = -jr/2 + 7$. Furthermore, 
the dimensionless parameters 11^5 and nl 9 arc defined by 
n^/n! 9 = n th /(n th + l) and H*? + K e « = 1. 

Quantum detailed balance condition: In order to check 
condition Eq. 1471) . we note that the evolutions Eq. 153|) 
and H54f> can be cast in the superoperator form 



L(u)[« 



1 



1 + 2n th 



K(u)£ th [.] + ^4^ C ^ ( 57 ) 



where K^u) = iiT<i>(u) 
tionary state reads 



K(u)/2. The corresponding sta- 



ps =n?|+> 



n e _ 9 |-)(-| 



(58) 



which due to the time reversal invariance of Hamiltonian 
eigenvectors satisfies p^ — p^ . Then, it is easy to prove 
that Eq. (|47|l is satisfied identically. Consistently, notice 
that the underlying Markovian dynamic Eq. I|49|) satisfies 
the conditions Eq. fZty . 

Quantum regression theorem: As the quantum detailed 
balance condition is satisfied, the QRT is valid in an 
asymptotic regime. Consistently the stationary state of 
Eq. l|4^|) does not depend on 7^ [p^ = p^]. 

The transient deviation from the QRT can be easily 
obtained for this example. First, we note that the density 
matrix evolution defined by Eq. <|57|l is equivalent to the 
non-Markovian Bloch equation 



dSx(t) 
dt 

dS Y (t) 
dt 

dSzjt) 
dt 



- [ dTK*(t-T)S x (T), (59a) 
Jo 

- [ drK 9 (t -t)S y (t), (59b) 
Jo 

- f drK(t-r)[S z (T)-S^}, (59c) 
Jo 



where Sj(t) = Trg{ps(t)<Tj} are the expectation values 
of the Pauli matrixes aj, and <S|P = — 111 9 . In order to 
deal with diagonal matrixes, we analyze the correlations 
in the base A = {a x , o~ y , (a z — iS|?),I}. Then, the propa- 
gator for operator expectation values, A(£) = G(i)A(0), 
can be written as 

<&(*) = diag{P*(t),P*(t),ifc(*),l}- (60) 

Here, we defined the functions Pn(u) = [u + K^u)]^ 1 and 
P$(u) = [u + -ft'$(u)]~ 1 , which in term of the random rate 
set can be written in the time domain as 

P n (t) = (exp[- 7Jl t]) , P*(t) = e~^ t P Yl {t/2). (61) 

On the other hand, the extra inhomogeneous term 
[Eq. 1221] that defines the operator correlations, 
0(t)A(t + t) = G(T)0(t)A(t)+F(t, r), can be written as 

F(t,r) =G n (t,r)F n + G$(i,r)F$, (62) 

where we have defined the vectors 

F n = Tv s [O(0)A(0){p+(Q)-pf}}, (63a) 
F* = Tr s [O(0)A(0)p s (0)], (63b) 

with p±(0) = [p s (0) ± a z p s {0)cr z }/2 H|. We note that 
Fn measure the departure of the initial populations from 
the equilibrium values 11^? , while F$ measure the depar- 
ture of the initial coherences from their null stationary 
value. Thus, F(t, r) vanishes if the system start in the 
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equilibrium state p^. On the other hand, the time de- 
pendence of F(£,t) is defined by the matrixes 

Gn(t,r) = diag{/ o (i,r),/ (i,T),/n(i,T),0}, (64a) 
&*(i,r) = dmg{U{t,T),U(t,r)J o (r,t),0}, (64b) 

with the definitions 

/ (i,r) = e ^ T Pn(t + T/2)-PnW^(r), (65a) 
/*(i,r) = P»(t + r)-P*(t)P*(r), (65b) 
/n(i,r) = Pn(t + r)-Pn(t)Pn(r). (65c) 

These functions measure the transient departure from the 
validity of the QRT. Only when the decay behaviors are 
exponential, they vanish identically and the QRT is valid 
at all times. This situation happens when the evolution 
is Markovian. 



0.6- 




0.0 -, ' ' ■ ' ■ ' 




B. Transient decay behaviors 

In order to illustrate the previous results, we spec- 
ify the properties of the complex environment, which in 
the context of the GBMA means to characterize the set 
{7r,Pr}- We choose 

(1 - er a ) 

7 fl = 70 exp[-bR], P R = _ e „ ajV ^ exp[-aP], 

(66) 

where R S [0, N — 1], 7 scale the random rates, and the 
dimensionless constants b and a measure the exponen- 
tial decay of the random rates and their corresponding 
weights. The relevant parameters of this set are 

\1R) b 

Here, 7 is the average rate and 8 measures the dispersion 
of the random rate set. On the other hand, in the limit 
N — > 00 the set Eq. I|66|l may leads to system dynamics 
characterized by a power law behavior whose exponent 
is given by a 36] . 

In Fig. we plot the transient decay behavior of the 
correlation 



C X Y(t,T)=<j x (t)a y (t + T), (68) 

which from Eq. 16UI) and (|62|l can be written as 

C XY (t, T ) = t{P*(T)S z (t) + f o (t,T){S z (0)-S?]}, (69) 

with S z (t) = Sf + P n (t) [S z (0) - Sf) . We have chosen a 
cero temperature reservoir, nth = 0, characterized by the 
random rate set Eq. (|66|l . As initial condition we take the 
pure state |+) . Thus, Sf = -1 and S z (0) = 1. Notice 
that the initial value of each plot describe the decay of 
the initial condition from the upper to the lower state. 
In fact C XY (t,Q) = iS z (t). 



FIG. 1: Transient decay behavior of C' XY (t, t) = CxY(t,r)/i. 
The parameters of the complex environment are b — 2.15, 
a = ab, a = 1/2, TV = 5, nth = 0, and 7^/7 = 0.02. The 
dispersion rate results /3/j = 0.4. The dotted lines correspond 
to the QRT. From top to bottom, we set jt = 0.25, 0.75, 2.5, 
and 250. 



As can be seen from the graphics, the predictions of the 
QRT are asymptotically valid in the stationary regime, 
where the function /o(i, t) vanish identically. In fact, the 
correlation behavior predicted by the QRT follows from 
Eq. ljtj9"|) after replacing / (t, r) — > 0. 

The transient deviations from the QRT are propor- 
tional to the departure of the system decay behavior from 
an exponential one. This departure arises from the com- 
petence between the exponential decay introduced by the 
rate 7$ and the non-Markovian effects induced by the 
random rate dispersion. From Eq. (|69|) it is evident that 
the dispersive rate 7$ introduces a global exponential 
decay. Thus, in general, by increasing this rate, the tran- 
sient deviation from the QRT are diminished. On the 
other hand, an increasing of (3 implies a strong deviation 
from an exponential decay. 

In order to enlighten the dependence in the dispersion 
of the random rate set, in Fig. [21 we plot P$(t), Eq. 16 H . 
for different values of the dispersion r ate 8 . This function 
determine both the coherence decay [46j and the devia- 
tions /o(t, t) and /$(t, r), Eq. (|d5|) . 

The short time behavior can be approximated by the 
exponential decay P$(t) ~ exp[— {7$ + (■yn)/2}i\, while 
the asymptotic one by P$(t) — exp[— {7$ + (l/jn)/2}t]. 
These behaviors can be straightforwardly obtained from 
Eq. I|61l) . In the intermediate regime the decay is approx- 
imately a power law with exponent a. By diminishing 8, 
the non-exponential decay behaviors occurs at small val- 
ues of P$(i). In fact, for 8/-f <C 1, the whole decay may 
be well approximated by P$(i) ~ exp[— {7$ + {"fR)/2}t] 
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1000 



FIG. 2: Decay behavior of P$(t). From top to bottom, the 
parameters of the complex environment are b = 2.15, 6.05, 
10.6 and 15.2. In all cases we take a — ab, a = 1/2, N = 5, 
nth — 0, and 7$/7 = 0.02. The dotted line corresponds to the 
Markovian decay exp[— i(7$+7/2)]. The inset shows the same 
curves in a linear plot. In this scale, the decay for b — 10.6 
and 15.2 are indistinguishable from the exponential one. 



(see inset). Thus, in this situation, the QRT may be as- 
sumed valid at all times for correlations involving the de- 
viations fo(t, r)«0 and r) ~ 0. On the other hand, 
as the population decay Pa(t) |4g does not involves the 
dispersive rate 7$, in general we can not disregard the 
transient effects introduced by fn(t,r), Eq. (|65cl) . 



C. Decay under the action of an external field 



For dealing with a manageable dynamics, we consider 
the external Hamiltonian H ext (t) = {Ml/2)(a^ e~ iuJAt + 
ae +iuj A ty Then, the system density matrix dynamics can 
be associated with a spin subject to a resonant external 
magnetic field [j| or with a two level optical transition 
driven by a resonant laser field Q. 

In an interaction representation with respect to 
hojA&z/2, the effective system Hamiltonian reads -ffg = 
hVLcr x /2. Thus, the evolution of the states pn{t) is given 
by Eq. |@SJ with H s — ► H e J f (see Appendix A). From 
Eq. H2bafl and (|27[) . the expectation values of the Pauli 
matrixes evolve as 



dS x (t) 
dt 

dS Y (t) 
dt 

dS z (t) 
dt 



drT x {t-r)S x (T), 



(70a) 



-nS Z (t)-J dT{T Y {t-T)S Y {T) 

+T(t-T)[Sz(r)-S^}}, (70b) 

ns Y {t) + J dT{r{t-T)s Y {r) 

-T z {t-T)[S z {r)-S%]}. (70c) 



In Appendix A we give the exact expressions for the ker- 
nels T j(t), j — x, y, z, and T(i), as well as the expression 
for the non local superoperator L(it), Eq. (|14(l . We re- 
mark that independently of the set of random rates {7_r} 
and weights {Pr}, the kernels that define the evolution 
Eq. (|70|l depend explicitly on the intensity parameter $1. 

The stationary state corresponding to the evolution of 
each state pn(t), Eq. (|4*9)l . reads 



n°° 



1 



I 



(71) 



2 f ' (l + 2n^)[ 7i?7 * + n 2 ] 
which explicitly depends on 7^ if Q, ^ 0. Then, when 



1.0 - 



0.0 



p/y=5x10 




y /y=0.02 £2/y=0.2 
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FIG. 3: Stationary decay behavior of C||(t,r). From top to 
bottom, the parameters of the complex environment are b = 
10.6, 6.05, and 2.15. In all cases we take a = ab, a — 1/2, 
N = 5, nth = 0, and 7^/7 = 0.02. The intensity is fi/7 = 0.2. 
The dotted lines correspond to the QRT. 

the system is subject to the action of the external field 
the QRT in not fulfilled, even in the asymptotic regime. 
Consistently, as was demonstrated in Ref. the un- 
derlying Markovian evolution Eq. (|49(l does not satisfy 
the detailed balance condition Eq. I|44|) . neither the su- 
peroperator L(u) [Eq. HA10|) ] satisfy Eq. g7|). As the 
QRT is not fulfilled, the operators correlations must to 
be calculated from the microscopic Hamiltonian dynam- 
ics, which in our case implies the averaging procedure 
corresponding to the GBMA. 

In the next figures we characterize the correlation 



C U (t,r) 



{C xx (t,T) + C YY (t,T)}/4 (72) 
-i{C XY (t,T)-C YX (t,r)}/4, 



where Cjk(t,r) = <jj(t)ak(t + r) are the correlations 
of the Pauli matrixes. Then it follows C^i(t,T) = 
a^(t)a(t + t). Each contribution Cjk(t,r) can be deter- 
mine from Eq. I|33fl . which involves an average of the 
corresponding Markovian solutions over the random rate 
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set, Eq. H66|) . On the other hand, the QRT predictions 
follows from Eq. ® with F(t, r) -> 0. 

In Fig. |2| we plot the stationary decay 
C u (oo,t)/C u (oo,0), where C u (oo,0) = [1 + 5 z (oo)]/2, 
for different values of the rate (3. We note that both 
the decay behaviors and stationary values differ from 
the QRT predictions. As the evolution of Sx(t) does 
not depend on ft [see Eq. I|7(}a|l ]. in the asymptotic 
regime Cxx(t,r) satisfies the QRT. Furthermore, as 
\im t ->oo Sx(t) = 0, from Eq. l(3l)l) and il3D|) we deduce 
that the disagreement in the asymptotic values with 
respect to the QRT only arises due to the contribu- 
tion Cyy(f, r), while CxY{t, r) and Cyx{t, r) only 
contribute to the difference in the decay behaviors. 




20 40 60 80 

y x 
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FIG. 4: Stationary decay behavior of C-\i(t, r) after increasing 
the dispersive rate, 74-/7 = 0.1. The remaining parameters 
are the same than in Fig. [3] The dotted lines correspond to 
the decay predicted by the QRT. 



The asymptotic value predicted from Eq. is 



C u (oo,oo) 



/ ^W2 \ 2 



(l + 2n th ) 2 
while from Eq. iffify. for the QRT we get 



C n {oo, 00) 



QRT 



(1 



2n t?l ) 2 ( 



7R7i 



tt 2 



(73) 



(74) 



where we have used the stationary state Eq. Q7ip. As 
can be seen in the graphics, the difference between both 
predictions grows by increasing the dispersion rate j3. 

In Fig. 2| , we plot the same correlation after increasing 
the dispersive rate 7$ and maintaining fixed all other 
parameters. We note that the deviations with respect to 
the QRT are diminished. In fact, for small values of /3/7, 
the dynamical deviations goes asymptotically to zero. 

In Fig. we plot C^i(oo,t) for different values of the 
field intensity fl. The deviations with respect to the QRT 



FIG. 5: Stationary decay behavior of Cn(t,r) for different 
values of the intensity parameter. From top to bottom, we 
take fi/7 = 1 and 5. In both cases the parameters of the 
complex environment are b = 2.15, a — ab, a — 1/2, TV = 5, 
nth — 0, and 74-/7 = 0.1. The dotted lines correspond to the 
decay predicted by the QRT. 



are diminished by increasing f2. Even more, in the limit 
of high intensity, the dynamical deviations vanish. 

The previous parameter dependence analysis relies in 
a specific correlation and random rate set model. Similar 
conclusion can be obtained by studying the asymptotic 
behaviors predicted by Eq. (|36|l and Eq. I|39|l for arbi- 
trary correlations and random rate sets. The deviations 
between both equations are proportional to the random 
dispersion of the stationary state p^, Eq. ifTTl) . Thus, as 
a measure of the departure from the validity of the QRT 
in the stationary regime, we introduce the matrix 



(75) 

(Ft) 

with j, k — x,y,z, and where (00) = Trs{p|jr tjj\ 
are the random stationary Pauli expectation values. A 
general characterization of this matrix can be given in a 
small and high intensity limits. 

In the small intensity limit, <C {tr}> w e can 
approximate 



n 2 



(1 + 2n th ) 2 



{<K) 2 )-<r|) 2 } + 0(^), (76) 



R 



(7fl/2 



7$) 1 = 1/7|. Consistently, E JK 



where r 

goes to zero in the limit of small intensity £1. On the 
other hand, by increasing the dispersive rate 7$, each 
contribution in Eq. I|7()|l diminish, which in turn means 
that the predictions of the QRT approach the exact dy- 
namics. 

In the high intensity limit, f2 3> {7i?} we get 



n- 



(f + 2n th ) 



;{(il)-{iR?} + o{n- 



(77) 
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This expression implies that by increasing fi, the validity 
of the QRT is asymptotically recuperated. This result is 
consistent with the fact that at high intensity values jljj 
the underlying Markovian dynamics for pn(t) satisfies the 
detailed balance condition Eq. (|44ll . In fact, in this limit 
the stationary states p^ [Eq. (|71J) ] can be approximated 
by ss 1/2, which as expected does not depends on the 
random rate. 



VI. SUMMARY AND CONCLUSIONS 



by a non-local Lindblad evolution. 
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In this paper we obtained the conditions under which a 
QRT can be assumed valid for quantum non-Markovian 
master equations defined by Lindblad superoperators 
with memory elements. In order to work on the base of a 
full Hamiltonian description we deduced our results from 
a GBMA. This approximation in a natural way leads to 
these kind of equations. In this context, we demonstrated 
that operator correlations follow from a weighted average 
of a set of Markovian solutions, each one characterized 
by a different dissipative rate. 

From our analysis, we deduced that a non-Markovian 
QRT can only be granted in an stationary regime if 
the evolution satisfies a non-Markovian quantum detailed 
balance condition [Eq. (|47|l ]. which in turn is related with 
the time reversal symmetry of the stationary dissipative 
dynamics. When this is not the case, the QRT is not 
fulfilled at any time, and in consequence, the only way 
of calculating operators correlations is from the corre- 
sponding microscopic dynamics. We remark that the im- 
possibility of formulating a non-Markovian quantum re- 
gression theorem outside a stationary regime can be also 
demonstrated from general dynamical arguments (see 
Appendix B). 

In general, the departure from the predictions of the 
QRT not only implies differences in the decay behaviors, 
but also in the asymptotic values of the operators corre- 
lations. The magnitude of these deviations are propor- 
tional to the departure of the system dynamics from a 
semigroup dynamical behavior, i.e., an exponential one. 

As an example of our results we worked out the dynam- 
ics of a two level system subject to the action of an ex- 
ternal coherent field and a complex thermal environment 
whose action can be described in a GBMA. Without the 
external field, the QRT is valid in an asymptotic regime. 
Consistently, the non-Markovian quantum detailed bal- 
ance condition is also satisfied. The presence of the ex- 
ternal field invalidates the QRT, even in the stationary 
regime. Nevertheless, in the limit of high intensity, or 
when the effect of a Markovian dispersive contribution 
is dominant, the QRT is asymptotically reestablished to 
lowest order in the corresponding expansion parameters. 

The present results provides an step forward in the 
understanding of non-Markovian open quantum systems 
dynamics. In fact, we have found solid physical criteria 
for the possibility of using a QRT for calculating opera- 
tors correlations when the system dynamics is described 



APPENDIX A: DENSITY MATRIX EVOLUTION 

Here we characterize the evolution of the density ma- 
trix elements for the example developed in Section V. 



1. Markovian evolution 



By denoting the matrix elements as 

( nf>(t) $f>(t) 



(Al) 



in an interaction representation with respect to Hloao' z /2, 
from Eq. I|49|) we get the evolutions 



l4 fl) (t) = lR {TU e _ q UX t >(t)±U^U w (t)} (A2a) 



dt 



±f[C«'i -<i>L'"./)]. 

d _ (m , « . if2 , 



f f ^\t) = ±f[n^>(t)-U w (t)}- 1 %^>(t).(A2b) 

The operators expectation values are defined by the evo- 
lutions 



r(«)m_Tr( fi )/ 



dS ( x R \t) 
dt 

,(.R), 



= -its^\t), 



(A3a) 



d -^l = -nsi R \t)- 7 *s¥\t), (A3b) 

= nS^(t)- lR [si R \t)-Sn (A3c) 

where Sj (t) = Tis{pR(t)&j}, with <jj the Pauli ma- 
trixes. 



2. Non-Markovian evolution 

Here, for arbitrary set {jr,Pr}, we present the ex- 
act expressions for the kernels that define the evolution 
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Eq. (25- From Eq. gJJ and Eq. l|A3)) . we get 



(A4a) 



r P 



IY(u) = Z?{( U + C)[^ + (u + 7<I ,)]+r! 2 } + 7 $ ,(A4b) 
r z («) = 2D{(u + B)[^ + {u + lq> )] + Q 2 }, (A4c) 



T(u) = D(C-B)Q, 

where D denotes the function 

B(u)/2 



D(u) 



(A4d) 



(A5) 



[u + B{u)][u + B(u)/2 + 7$] + tt 2 ' 
The extra function B and C are defined by 



<n«)> ' 

where we have introduced 
T(u) = 



{T(u)^r} 



1/2 



{u + lR )[u + lR /2 + ^] + n 2 ' 



(A7) 



The corresponding density matrix evolution can be writ- 
ten as 



ri±N = r z ( u ){Tni 9 n + ( u )±n + 3 n_( u )}, 

iCl(u) . 



±- 



[*+(«)-*-(«)]. 



(A8a) 



± / \ r „ , . „ . T(u)/(2u) 

$ ± u = ±-^pi+(« -n_ (u )}±i y ' /K ' 

2 1 + 2n th 



-r+(tt)$±(u)-r»(u)$ T («), 



(A8b) 



where, for shortening the notation, = ug(u) — g(0) 
denotes the Laplace transform of the time derivative of 
an arbitrary function g(t). In the inhomogeneous term 
for the coherences, we have used 11+ (it) + U-(u) = 1/u. 
Furthermore, we have defined Cl(u) = Q + T(u) and 
= ^[Tx(u) ± Ty(u)]. The stationary state reads 

_ i f nr z a y - [T Y T Z + T(T + n)]a z \ 

Ps ~2\ l+ (i + 2n t/l )[r y r z + (T + ^] /' (A9) 

with the notation Tj = Fj(u)\ u=0 . Consistently, after 
some algebra, it is possible to write this state as an av- 
erage of the corresponding Markovian stationary states, 
i.e., pf = (pf ), where p^ is defined by Eq. (J7TJ- 

The superoperator L(u) [Eq. I|14|l ] corresponding to 
the evolution Eq. (|A8|> can be written as non-diagonal 
non-local Lindblad superoperator 

L(«)H = C H (u)[.] + a "P(u)([V a , *Vl] + [V a ; V$]), 

a/3 

(A10) 

with the operators {y a } Q =i, 2, 3 = {o~, a\ & z }. The Hamil- 
tonian contribution reads 



C H (u)[*]=-i^P-[a x ,.} 



(All) 



and the matrix elements a a p(u) are defined by 



033 (" 

012 (U 

a 13 (u 

031 (U 



r Y (u) 



a 2 i{u) = --{T x (u 



rz(«)}, 
-!>(«)}, 



a 2 3(w) 
a 32 (w) 



'4(l + 2n th )' 
T(«) 



(A12) 
(A13) 

(A14) 
(A15) 
(A16) 

(A17) 



4(l + 2ntfc) 

Without the external excitation, 17 = 0, the super- 
operator L(it) reduce to Eq. I|57[l. Furthermore, when 
the coherence decay behavior can be approximated by 
an exponential one, P$(i) = e~ 1 * t Pn(t/2) — exp[— (7$ + 
("fn)/2)t\, the density matrix evolution can be written in 
a Schrodinger representation as 

d Ps (t) 



dt 



±[H s , Ps (t)] + ^U[p s (t)) (A18) 



+ 



1 



1 + 2n th 



f dTK(t-r)C th [ps(r)}, 
Jo 



with Hs = Huja<J z /2. This expression relies in the valid- 
ity of the approximation K(u ± iwa) — K(oo) = {jr), 
which can be considered always valid if wa is an opti- 
cal frequency. Furthermore, if 7$ <c (jr) the dispersive 
contribution can be drop. In general this last condition 
is valid when the decay of Pn(t) develops two strong dif- 
ferent time scales. For example, consider a random rate 
that assumes only two different values Jyi, with prob- 
abilities P T/i . Then P n {t) = P T e~ 7 T« + p^'H*. Under 



the conditions Pi -C Pf 
approximate Pri^e^ 7 ** ; 



and j 1 <C 7$ <C 7 f, we can 
i e -< 7 «>* + 0(P x /P T ). Another 
examples follow from the decay of Fig. [3 for small /3/ 7 . 
On the other hand, Eq. I|A18|) can also be assumed valid 
in presence of the external field if the exact kernels are 
taken to cero order in the intensity parameter £1. 



APPENDIX B: ON THE IMPOSSIBILITY OF 
FORMULATING A NON-MARKOVIAN 
QUANTUM REGRESSION THEOREM AT ALL 
TIMES 

The impossibility of formulating a non-Markovian re- 
gression theorem outside a stationary regime can be 
demonstrated on general dynamical arguments. In fact, 
it is simple to proof that the validity of the quantum re- 
gression theorem at all times is only compatible with a 
Markovian dynamics. This affirmation seems to contra- 
dict our main conclusions. Nevertheless, here we demon- 
strate that this result confirm the correctness of our ap- 
proach. 

First, we write the system density matrix as 



p s (t) = T(t)[p s (0)}, 



(Bl) 
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where T(t) is the propagator corresponding to the evolu- 
tion Eq. H13|) . Then, it is defined in the Laplace domain 
by 



T(u) 



- [C H + L(u)] ' 



(B2) 



In terms of this object, we can write the operator expec- 
tation values as 



A(t) = Tr s {A(0)T(t)[ (Os (0)]}, (B3a) 
= Tr s { Ps (0)T#(t)[A(0)]}, (B3b) 

where the second line defines the dual propag ator T*(t). 
By assuming valid the quantum regression theorem, the 
operator correlations can be written as 0, IE IS 13 

0(t)A(t + r) = Tr s { Ps (t)O(0)T#(r)[A(0)]}. (B4) 

This expression must to be valid for arbitrary operators 
O and A. In particular, by taking O = Is, where Is is 
the system identity operator, it follows 



A(t + r) = Tr s {ps(t)T*(T)[A(0)}}, (B5a) 
= Tr s {A(0)T(r)[p s (i)]}, (B5b) 
= Trg{A(0)T(T)T(t)[p g (0)]}. (B5c) 



Eq. IjBlfl . this equality can only be satisfied if the prop- 
agator T(t) corresponds to a semigroup structure, i.e., 
a Markovian evolution. Therefore, a regression theorem 
can be satisfied at all times only when the dynamics does 
not has any memory contribution. We notice that this 
result is in perfect agreement with our main conclusions. 
In fact, we have found that a non-Markovian quantum re- 
gression theorem may be valid (or not) only in a station- 
ary regime. In this limit, the previous calculations steps 
does not impose any constraint on the propagator T(i). 
This affirmation follows trivially by taking ps(0) = Pg in 
Eq. (|B7|I . or equivalently by introducing the limit t — > oo, 



HmT(t + t) Ps (0) = hm T(r)T(t)p s (0), (B8) 



which, independently of the properties of T(t), deliver 
Ps° ~ Ps°- This last equality follows immediately from 
T(r) [/j^P] = pg , expression valid for any time r. Alterna- 
tively, one can take the limit t — > oo in Eq. (|B5c|l 



lim A(t + t) 

t — >oo 



On the other hand, from Eq. l|B3a(l . we can write 



Tr 5 {A(0)T(r) lim T(t)[p s (0)]}, (B9a) 

t—*oo 

= Tr s {A(0)T(r)[^]}, (B9b) 

= Tr s {A(0)^}. (B9c) 
On the other hand, in the same limit, from Eq. (|B6I) . as 

expected, we get the same result 



A(t + r) = Tr s {A(0)T(i + r)[p s (Q)}}. (B6) 

As A(0) is an arbitrary operator, by comparing this ex- 
pression and Eq. (|B5c|l . it follows 



(B7) 



T(t + t) Ps (0) =T(T)T(t) Ps (0). 



For arbitrary time t < oo, and ps(0) ^ p]?, where pg 
is the stationary state corresponding to the dynamics 



lim A(t + r)=Tr s {A(0)^}. 



(BIO) 



Therefore, the calculations steps that lead to the con- 
straint Eq. (|B7(1 only contradict the possibility of es- 
tablishing a non-Markovian quantum regression theorem 
outside the stationary regime. These arguments provide 
an alternative demonstration of the consistency and cor- 
rectness of our results. 
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